Model11test
Several authors
Created: 12-04-2024. Last modified: 30-04-2024.
Let us set some global options for all code chunks in this document.
knitr::opts_chunk$set(
message = FALSE, # Disable messages printed by R code chunks
warning = FALSE, # Disable warnings printed by R code chunks
echo = TRUE, # Show R code within code chunks in output
include = TRUE, # Include both R code and its results in output
eval = TRUE, # Evaluate R code chunks
cache = FALSE, # Enable caching of R code chunks for faster rendering
fig.align = "center",
out.width = "100%",
retina = 2,
error = TRUE,
collapse = TRUE
)
rm(list = ls())
set.seed(1982)1 Preprocessing
Let us now load some required libraries.
# Load required libraries
# inla.upgrade(testing = TRUE)
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)
inla.setOption(num.threads = 7)
library(inlabru)
library(rSPDE)
library(MetricGraph)
library(plotly)
library(dplyr)
library(sf)
library(here)Function standarize() below is later used to standardize
the covariate SpeedLimit.
standardize <- function(x) {return((x - mean(x)) / sd(x))}To keep track of the changes, we provide summaries of every new created object. Those summaries can be accessed by pressing the Show buttons below
We load the graph object sf_graph (which only contains
weights) and the data (already graph-processed).
load(here("Graph_objects/graph_construction_25_04_2024partialtomtomwhichlonglatsf.RData"))
load(here("Data_files/data_day7142128_hour13_with_no_consecutive_zeros_partialtomtom_graph_25_04_2024_processed.RData"))
data_on_graph = data_on_graph %>%
dplyr::select(-datetime)sf_graph$get_edge_lengths() %>% head() %>% capture.output() %>% grep("^Units:", ., value = TRUE)
## [1] "Units: [km]"summary(sf_graph)
## A metric graph object with:
##
## Vertices:
## Total: 4017
## Degree 2: 1821; Degree 3: 180; Degree 4: 1402; Degree 5: 94; Degree 6: 367;
## Degree 7: 36; Degree 8: 116; Degree 12: 1;
## With incompatible directions: 0
##
## Edges:
## Total: 6827
## Lengths:
## Min: 0.002834658 ; Max: 0.2743201 ; Total: 311.1441
## Weights:
## Columns: Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour
## That are circles: 0
##
## Graph units:
## Vertices unit: degrees ; Lengths unit: km
##
## Longitude and Latitude coordinates: TRUE
## Which spatial package: sf
## CRS: EPSG:4326
##
## Some characteristics of the graph:
## Connected: TRUE
## Has loops: FALSE
## Has multiple edges: TRUE
## Is a tree: FALSE
## Distance consistent: TRUE
## Has Euclidean edges: FALSE
##
## Computed quantities inside the graph:
## Laplacian: FALSE ; Geodesic distances: TRUE
## Resistance distances: FALSE ; Finite element matrices: FALSE
##
## Mesh: The graph has no mesh!
##
## Data: The graph has no data!
##
## Tolerances:
## vertex-vertex: 0.001
## vertex-edge: 0.001
## edge-edge: 0
summary(data_on_graph)
## ID speed day .distance_to_graph
## Min. :5701 Min. : 0.000 Min. :1.000 Min. :0.000000
## 1st Qu.:6571 1st Qu.: 1.609 1st Qu.:2.000 1st Qu.:0.001655
## Median :6687 Median :14.484 Median :3.000 Median :0.003552
## Mean :7106 Mean :15.142 Mean :2.556 Mean :0.004469
## 3rd Qu.:7281 3rd Qu.:24.140 3rd Qu.:4.000 3rd Qu.:0.006152
## Max. :8969 Max. :99.779 Max. :4.000 Max. :0.019991
## .edge_number .distance_on_edge .group .coord_x
## Min. : 1 Min. :0.0000 Length:39574 Min. :-122.4
## 1st Qu.:1413 1st Qu.:0.2726 Class :character 1st Qu.:-122.4
## Median :2879 Median :0.5149 Mode :character Median :-122.4
## Mean :3090 Mean :0.5089 Mean :-122.4
## 3rd Qu.:4753 3rd Qu.:0.7501 3rd Qu.:-122.4
## Max. :6826 Max. :1.0000 Max. :-122.4
## .coord_y
## Min. :37.77
## 1st Qu.:37.78
## Median :37.79
## Mean :37.79
## 3rd Qu.:37.79
## Max. :37.81The following commands remove zero speed observations that are 1m away from the graph, and after that, they remove any speed observations that are 3m away from the graph.
to_remove = data_on_graph %>%
filter(speed == 0, .distance_to_graph > 0.001)
data_on_graph = setdiff(data_on_graph, to_remove) %>%
filter(.distance_to_graph <= 0.003)summary(to_remove)
## ID speed day .distance_to_graph .edge_number
## Min. :5701 Min. :0 Min. :1.000 Min. :0.001001 Min. : 1
## 1st Qu.:6574 1st Qu.:0 1st Qu.:2.000 1st Qu.:0.002757 1st Qu.:1376
## Median :6688 Median :0 Median :3.000 Median :0.004563 Median :2788
## Mean :7078 Mean :0 Mean :2.556 Mean :0.005453 Mean :3031
## 3rd Qu.:7277 3rd Qu.:0 3rd Qu.:4.000 3rd Qu.:0.007129 3rd Qu.:4715
## Max. :8969 Max. :0 Max. :4.000 Max. :0.019988 Max. :6812
## .distance_on_edge .group .coord_x .coord_y
## Min. :0.0000 Length:8597 Min. :-122.4 Min. :37.77
## 1st Qu.:0.3078 Class :character 1st Qu.:-122.4 1st Qu.:37.78
## Median :0.5301 Mode :character Median :-122.4 Median :37.79
## Mean :0.5177 Mean :-122.4 Mean :37.79
## 3rd Qu.:0.7362 3rd Qu.:-122.4 3rd Qu.:37.79
## Max. :1.0000 Max. :-122.4 Max. :37.81
summary(data_on_graph)
## ID speed day .distance_to_graph
## Min. :5701 Min. : 0.000 Min. :1.00 Min. :0.0000000
## 1st Qu.:6565 1st Qu.: 9.656 1st Qu.:2.00 1st Qu.:0.0005835
## Median :6683 Median :19.312 Median :3.00 Median :0.0012618
## Mean :7115 Mean :19.385 Mean :2.55 Mean :0.0013511
## 3rd Qu.:7286 3rd Qu.:27.359 3rd Qu.:4.00 3rd Qu.:0.0021077
## Max. :8969 Max. :99.779 Max. :4.00 Max. :0.0029999
## .edge_number .distance_on_edge .group .coord_x
## Min. : 1 Min. :0.0000 Length:14534 Min. :-122.4
## 1st Qu.:1376 1st Qu.:0.2552 Class :character 1st Qu.:-122.4
## Median :2956 Median :0.5139 Mode :character Median :-122.4
## Mean :3130 Mean :0.5099 Mean :-122.4
## 3rd Qu.:4803 3rd Qu.:0.7697 3rd Qu.:-122.4
## Max. :6817 Max. :1.0000 Max. :-122.4
## .coord_y
## Min. :37.77
## 1st Qu.:37.78
## Median :37.79
## Mean :37.79
## 3rd Qu.:37.79
## Max. :37.81We add data to the graph.
sf_graph$add_observations(data = data_on_graph,
group = "day",
normalized = TRUE,
clear_obs = TRUE)sf_graph$get_data()
## # A tibble: 14,534 × 9
## ID speed day .distance_to_graph .coord_x .coord_y .edge_number
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 6666 0 1 0.00100 -122. 37.8 1
## 2 8941 12.9 1 0.00230 -122. 37.8 4
## 3 8768 24.1 1 0.00233 -122. 37.8 6
## 4 8929 32.2 1 0.00151 -122. 37.8 6
## 5 8965 0 1 0.000647 -122. 37.8 9
## 6 8965 19.3 1 0.00170 -122. 37.8 14
## 7 8954 22.5 1 0.00103 -122. 37.8 14
## 8 8774 19.3 1 0.00279 -122. 37.8 14
## 9 6655 30.6 1 0.00208 -122. 37.8 18
## 10 6677 14.5 1 0.000436 -122. 37.8 20
## # ℹ 14,524 more rows
## # ℹ 2 more variables: .distance_on_edge <dbl>, .group <chr>
summary(sf_graph)
## A metric graph object with:
##
## Vertices:
## Total: 4017
## Degree 2: 1821; Degree 3: 180; Degree 4: 1402; Degree 5: 94; Degree 6: 367;
## Degree 7: 36; Degree 8: 116; Degree 12: 1;
## With incompatible directions: 0
##
## Edges:
## Total: 6827
## Lengths:
## Min: 0.002834658 ; Max: 0.2743201 ; Total: 311.1441
## Weights:
## Columns: Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour
## That are circles: 0
##
## Graph units:
## Vertices unit: degrees ; Lengths unit: km
##
## Longitude and Latitude coordinates: TRUE
## Which spatial package: sf
## CRS: EPSG:4326
##
## Some characteristics of the graph:
## Connected: TRUE
## Has loops: FALSE
## Has multiple edges: TRUE
## Is a tree: FALSE
## Distance consistent: TRUE
## Has Euclidean edges: FALSE
##
## Computed quantities inside the graph:
## Laplacian: FALSE ; Geodesic distances: TRUE
## Resistance distances: FALSE ; Finite element matrices: FALSE
##
## Mesh: The graph has no mesh!
##
## Data:
## Columns: ID speed day
## Groups: .group
##
## Tolerances:
## vertex-vertex: 0.001
## vertex-edge: 0.001
## edge-edge: 0We get the values of the weights at data locations. This essentially gives us covariates from the weights.
sf_graph$edgeweight_to_data(data_loc = TRUE)sf_graph$get_data()
## # A tibble: 57,108 × 54
## ID speed day .distance_to_graph Length FRC SpeedLimit StreetName
## <int> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 NA NA NA NA 0.0361 5 40 Harrison St
## 2 NA NA NA NA 0.0361 5 40 Harrison St
## 3 NA NA NA NA 0.0361 5 40 Harrison St
## 4 6666 0 1 0.00100 0.0361 5 40 Harrison St
## 5 NA NA NA NA 0.0361 5 40 Harrison St
## 6 NA NA NA NA 0.0361 5 40 Harrison St
## 7 NA NA NA NA 0.0361 5 40 Harrison St
## 8 NA NA NA NA 0.0361 5 40 Harrison St
## 9 NA NA NA NA 0.0361 5 40 Harrison St
## 10 8941 12.9 1 0.00230 0.112 6 35 Rhode Island St
## # ℹ 57,098 more rows
## # ℹ 46 more variables: harmonicAverageSpeed <dbl>, medianSpeed <dbl>,
## # averageSpeed <dbl>, sampleSize <int>, averageTravelTime <dbl>,
## # medianTravelTime <dbl>, travelTimeRatio <dbl>, List_Number <int>,
## # `5percentile` <int>, `10percentile` <int>, `15percentile` <int>,
## # `20percentile` <int>, `25percentile` <int>, `30percentile` <int>,
## # `35percentile` <int>, `40percentile` <int>, `45percentile` <int>, …
summary(sf_graph)
## A metric graph object with:
##
## Vertices:
## Total: 4017
## Degree 2: 1821; Degree 3: 180; Degree 4: 1402; Degree 5: 94; Degree 6: 367;
## Degree 7: 36; Degree 8: 116; Degree 12: 1;
## With incompatible directions: 0
##
## Edges:
## Total: 6827
## Lengths:
## Min: 0.002834658 ; Max: 0.2743201 ; Total: 311.1441
## Weights:
## Columns: Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour
## That are circles: 0
##
## Graph units:
## Vertices unit: degrees ; Lengths unit: km
##
## Longitude and Latitude coordinates: TRUE
## Which spatial package: sf
## CRS: EPSG:4326
##
## Some characteristics of the graph:
## Connected: TRUE
## Has loops: FALSE
## Has multiple edges: TRUE
## Is a tree: FALSE
## Distance consistent: TRUE
## Has Euclidean edges: FALSE
##
## Computed quantities inside the graph:
## Laplacian: FALSE ; Geodesic distances: TRUE
## Resistance distances: FALSE ; Finite element matrices: FALSE
##
## Mesh: The graph has no mesh!
##
## Data:
## Columns: ID speed day Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour
## Groups: .group
##
## Tolerances:
## vertex-vertex: 0.001
## vertex-edge: 0.001
## edge-edge: 0When running
sf_graph$edgeweight_to_data(data_loc = TRUE), some
NA values are created (because the data is grouped). We
remove them below. We also standardize the SpeedLimit
covariate.
data = sf_graph$get_data() %>%
drop_na(-StreetName) %>% # this drops all rows with at least one NA value but without taking into account StreetName
mutate(across(c("SpeedLimit"), ~standardize(.))) %>%
dplyr::select(speed, SpeedLimit)The code of chunk below was executed only one time.
{r, eval = FALSE}
aux = data |>
rename(distance_on_edge = .distance_on_edge, edge_number = .edge_number) |>
as.data.frame() |>
dplyr::select(edge_number, distance_on_edge, .group)
distmatrixlist = list()
for (i in 1:4) {
distmatrixlist[[i]] = sf_graph$compute_geodist_PtE(PtE = aux %>%
filter(.group == as.character(i)) %>%
dplyr::select(-.group),
normalized = TRUE,
include_vertices = FALSE)
}
save(distmatrixlist, file = here("Models_output/distmatrixfixed.RData"))The code of chunk above was executed only one time.
summary(data)
## speed SpeedLimit .group .edge_number
## Min. : 0.000 Min. :-2.3743 Length:14534 Min. : 1
## 1st Qu.: 9.656 1st Qu.:-0.1006 Class :character 1st Qu.:1376
## Median :19.312 Median :-0.1006 Mode :character Median :2956
## Mean :19.385 Mean : 0.0000 Mean :3130
## 3rd Qu.:27.359 3rd Qu.:-0.1006 3rd Qu.:4803
## Max. :99.779 Max. : 6.6173 Max. :6817
## .distance_on_edge .coord_x .coord_y
## Min. :0.0000 Min. :-122.4 Min. :37.77
## 1st Qu.:0.2552 1st Qu.:-122.4 1st Qu.:37.78
## Median :0.5139 Median :-122.4 Median :37.79
## Mean :0.5099 Mean :-122.4 Mean :37.79
## 3rd Qu.:0.7697 3rd Qu.:-122.4 3rd Qu.:37.79
## Max. :1.0000 Max. :-122.4 Max. :37.81We add the data again but now with the new standardized
SpeedLimit covariate.
sf_graph$add_observations(data = data,
group = "day",
normalized = TRUE,
clear_obs = TRUE)sf_graph$get_data()
## # A tibble: 14,534 × 7
## speed SpeedLimit .coord_x .coord_y .edge_number .distance_on_edge .group
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 0 -0.101 -122. 37.8 1 0.437 1
## 2 12.9 -0.617 -122. 37.8 4 0.144 1
## 3 24.1 -0.617 -122. 37.8 6 0.252 1
## 4 32.2 -0.617 -122. 37.8 6 0.658 1
## 5 0 -0.927 -122. 37.8 9 0.601 1
## 6 19.3 -0.927 -122. 37.8 14 0.0247 1
## 7 22.5 -0.927 -122. 37.8 14 0.362 1
## 8 19.3 -0.927 -122. 37.8 14 0.832 1
## 9 30.6 -0.927 -122. 37.8 18 0.358 1
## 10 14.5 -0.927 -122. 37.8 20 0.309 1
## # ℹ 14,524 more rows
summary(sf_graph)
## A metric graph object with:
##
## Vertices:
## Total: 4017
## Degree 2: 1821; Degree 3: 180; Degree 4: 1402; Degree 5: 94; Degree 6: 367;
## Degree 7: 36; Degree 8: 116; Degree 12: 1;
## With incompatible directions: 0
##
## Edges:
## Total: 6827
## Lengths:
## Min: 0.002834658 ; Max: 0.2743201 ; Total: 311.1441
## Weights:
## Columns: Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour
## That are circles: 0
##
## Graph units:
## Vertices unit: degrees ; Lengths unit: km
##
## Longitude and Latitude coordinates: TRUE
## Which spatial package: sf
## CRS: EPSG:4326
##
## Some characteristics of the graph:
## Connected: TRUE
## Has loops: FALSE
## Has multiple edges: TRUE
## Is a tree: FALSE
## Distance consistent: TRUE
## Has Euclidean edges: FALSE
##
## Computed quantities inside the graph:
## Laplacian: FALSE ; Geodesic distances: TRUE
## Resistance distances: FALSE ; Finite element matrices: FALSE
##
## Mesh: The graph has no mesh!
##
## Data:
## Columns: speed SpeedLimit
## Groups: .group
##
## Tolerances:
## vertex-vertex: 0.001
## vertex-edge: 0.001
## edge-edge: 0We build a mesh.
h = 0.05
sf_graph$build_mesh(h = h)summary(sf_graph)
## A metric graph object with:
##
## Vertices:
## Total: 4017
## Degree 2: 1821; Degree 3: 180; Degree 4: 1402; Degree 5: 94; Degree 6: 367;
## Degree 7: 36; Degree 8: 116; Degree 12: 1;
## With incompatible directions: 0
##
## Edges:
## Total: 6827
## Lengths:
## Min: 0.002834658 ; Max: 0.2743201 ; Total: 311.1441
## Weights:
## Columns: Length FRC SpeedLimit StreetName harmonicAverageSpeed medianSpeed averageSpeed sampleSize averageTravelTime medianTravelTime travelTimeRatio List_Number 5percentile 10percentile 15percentile 20percentile 25percentile 30percentile 35percentile 40percentile 45percentile 50percentile 55percentile 60percentile 65percentile 70percentile 75percentile 80percentile 85percentile 90percentile 95percentile road_type class_6 class_5 class_4 class_0 class_3 class_1 upto1 upto3 upto4 upto5 upto6 density density_per_hour
## That are circles: 0
##
## Graph units:
## Vertices unit: degrees ; Lengths unit: km
##
## Longitude and Latitude coordinates: TRUE
## Which spatial package: sf
## CRS: EPSG:4326
##
## Some characteristics of the graph:
## Connected: TRUE
## Has loops: FALSE
## Has multiple edges: TRUE
## Is a tree: FALSE
## Distance consistent: TRUE
## Has Euclidean edges: FALSE
##
## Computed quantities inside the graph:
## Laplacian: FALSE ; Geodesic distances: TRUE
## Resistance distances: FALSE ; Finite element matrices: FALSE
##
## Mesh:
## Max h_e: 0.04999798 ; Min n_e: 0
##
## Data:
## Columns: speed SpeedLimit
## Groups: .group
##
## Tolerances:
## vertex-vertex: 0.001
## vertex-edge: 0.001
## edge-edge: 0We get the value of the weights at mesh locations. This will allow us
to built matrices B.sigma and B.range below.
Again,
sf_graph$edgeweight_to_data(mesh = TRUE, add = FALSE, return = TRUE)
creates repeated information (because the data is grouped). We fix that
by filtering one group. We also standardize the SpeedLimit
covariate.
mesh = sf_graph$edgeweight_to_data(mesh = TRUE,
add = FALSE,
return = TRUE) %>%
filter(.group == 1) %>%
mutate(across(c("SpeedLimit"), ~standardize(.))) %>%
dplyr:::select.data.frame(SpeedLimit)summary(mesh)
## SpeedLimit
## Min. :-1.9800
## 1st Qu.:-0.1744
## Median :-0.1744
## Mean : 0.0000
## 3rd Qu.:-0.1744
## Max. : 5.16041.1 Stationary model
- Observe that we are considering replicates.
stat.time.ini <- Sys.time()
################################################################################
################################# STATIONARY MODEL #############################
################################################################################
rspde_model_stat <- rspde.metric_graph(sf_graph,
parameterization = "matern",
nu.upper.bound = 1.5)str(rspde_model_stat)
## List of 20
## $ f :List of 3
## ..$ model : chr "cgeneric"
## ..$ n : int 21507
## ..$ cgeneric:List of 5
## .. ..$ model: chr "inla_cgeneric_rspde_stat_general_model"
## .. ..$ shlib: chr "/usr/local/lib/R/site-library/INLA/bin/linux/64bit/external/rSPDE/librSPDE.so"
## .. ..$ n : int 21507
## .. ..$ debug: logi FALSE
## .. ..$ data :List of 5
## .. .. ..$ ints :List of 5
## .. .. .. ..$ n : int 21507
## .. .. .. ..$ debug : int 0
## .. .. .. ..$ graph_opt_i: int [1:135512] 0 0 0 0 0 0 0 0 0 0 ...
## .. .. .. ..$ graph_opt_j: int [1:135512] 0 1 230 984 985 986 1700 1701 2280 2455 ...
## .. .. .. ..$ rspde.order: int 2
## .. .. ..$ doubles :List of 11
## .. .. .. ..$ d : num 1
## .. .. .. ..$ nu.upper.bound : num 1.5
## .. .. .. ..$ matrices_less : num [1:94764] 0.0797 0 0 0 0 ...
## .. .. .. ..$ matrices_full : num [1:207848] 0.0797 0 0 0 0 ...
## .. .. .. ..$ start.theta : num [1:2] 0 0.223
## .. .. .. ..$ theta.prior.mean : num [1:2] 0 0.223
## .. .. .. ..$ prior.nu.loglocation: num -0.288
## .. .. .. ..$ prior.nu.mean : num 0.75
## .. .. .. ..$ prior.nu.prec : num 3
## .. .. .. ..$ prior.nu.logscale : num 1
## .. .. .. ..$ start.nu : num 0.75
## .. .. ..$ characters:List of 5
## .. .. .. ..$ model : chr "inla_cgeneric_rspde_stat_general_model"
## .. .. .. ..$ shlib : chr "/usr/local/lib/R/site-library/INLA/bin/linux/64bit/external/rSPDE/librSPDE.so"
## .. .. .. ..$ prior.nu.dist : chr "lognormal"
## .. .. .. ..$ parameterization : chr "matern"
## .. .. .. ..$ prior.theta.param: chr "theta"
## .. .. ..$ matrices :List of 2
## .. .. .. ..$ rational_table : num [1:5996] 999 6 0.001 0.412 0.005 ...
## .. .. .. ..$ theta.prior.prec: num [1:6] 2 2 0.1 0 0 0.1
## .. .. ..$ smatrices : list()
## .. ..- attr(*, "class")= chr "inla.cgeneric"
## $ cgeneric_type : chr "general"
## $ theta.prior.mean : num [1:2] 0 0.223
## $ prior.nu :List of 4
## ..$ loglocation: num -0.288
## ..$ mean : num 0.75
## ..$ prec : num 3
## ..$ logscale : num 1
## $ theta.prior.prec : num [1:2, 1:2] 0.1 0 0 0.1
## $ start.nu : num 0.75
## $ integer.nu : logi FALSE
## $ start.theta : num [1:2] 0 0.223
## $ stationary : logi TRUE
## $ rspde.order : num 2
## $ dim : num 1
## $ est_nu : logi TRUE
## $ nu.upper.bound : num 1.5
## $ prior.nu.dist : chr "lognormal"
## $ debug : logi FALSE
## $ type.rational.approx: chr "chebfun"
## $ mesh :Classes 'metric_graph', 'R6' <metric_graph>
## Public:
## add_mesh_observations: function (data = NULL, group = NULL)
## add_observations: function (data = NULL, edge_number = "edge_number", distance_on_edge = "distance_on_edge",
## build_mesh: function (h = NULL, n = NULL, continuous = TRUE, continuous.outs = FALSE,
## buildC: function (alpha = 2, edge_constraint = FALSE)
## buildDirectionalConstraints: function (alpha = 1)
## C: NULL
## characteristics: list
## check_distance_consistency: function ()
## check_euclidean: function ()
## clear_observations: function ()
## clone: function (deep = FALSE)
## CoB: NULL
## compute_characteristics: function (check_euclidean = FALSE)
## compute_fem: function (petrov = FALSE)
## compute_geodist: function (full = FALSE, obs = TRUE, group = NULL, verbose = 0)
## compute_geodist_mesh: function ()
## compute_geodist_PtE: function (PtE, normalized = TRUE, include_vertices = TRUE, verbose = 0)
## compute_laplacian: function (full = FALSE, obs = TRUE, group = NULL, verbose = 0)
## compute_PtE_edges: function ()
## compute_resdist: function (full = FALSE, obs = TRUE, group = NULL, check_euclidean = FALSE,
## compute_resdist_mesh: function ()
## compute_resdist_PtE: function (PtE, normalized = TRUE, include_vertices = FALSE, check_euclidean = FALSE,
## coordinates: function (PtE = NULL, XY = NULL, normalized = TRUE)
## drop_na: function (...)
## E: 1 3 3 6 4 7 9 11 13 14 15 16 18 19 20 21 23 25 27 29 31 ...
## edge_lengths: 0.0363234139144278 0.01586106867077 0.027923568765887 0. ...
## edges: metric_graph_edges
## edgeweight_to_data: function (loc = NULL, mesh = FALSE, data_loc = FALSE, weight_col = NULL,
## fem_basis: function (PtE)
## filter: function (..., .drop_na = FALSE, .drop_all_na = TRUE)
## geo_dist: list
## get_data: function (group = NULL, tibble = TRUE, drop_na = FALSE, drop_all_na = TRUE)
## get_degrees: function (which = "degree")
## get_edge_lengths: function (unit = NULL)
## get_edge_weights: function (data.frame = FALSE, tibble = TRUE)
## get_groups: function (get_cols = FALSE)
## get_initial_graph: function ()
## get_locations: function ()
## get_mesh_locations: function (bru = FALSE, loc = NULL, normalized = TRUE)
## get_PtE: function ()
## get_vertices_incomp_dir: function ()
## initialize: function (edges = NULL, V = NULL, E = NULL, vertex_unit = NULL,
## is_tree: function ()
## Laplacian: NULL
## mesh: list
## mesh_A: function (PtE)
## mutate: function (..., .drop_na = FALSE, .drop_all_na = TRUE)
## nE: 6827
## nV: 4017
## observation_to_vertex: function (tolerance = 1e-15, mesh_warning = TRUE)
## plot: function (data = NULL, newdata = NULL, group = 1, plotly = FALSE,
## plot_connections: function ()
## plot_function: function (data = NULL, newdata = NULL, group = 1, X = NULL, plotly = FALSE,
## plot_movie: function (X, plotly = TRUE, vertex_size = 5, vertex_color = "black",
## print: function ()
## process_data: function (data = NULL, edge_number = "edge_number", distance_on_edge = "distance_on_edge",
## prune_vertices: function (check_weights = TRUE, verbose = FALSE)
## PtV: NULL
## res_dist: NULL
## select: function (..., .drop_na = FALSE, .drop_all_na = TRUE)
## set_edge_weights: function (weights = rep(1, self$nE), kirchhoff_weights = NULL)
## summarise: function (..., .include_graph_groups = FALSE, .groups = NULL,
## summary: function (messages = FALSE, compute_characteristics = TRUE, check_euclidean = TRUE,
## V: -122.41277 -122.41249 -122.40376 -122.40358 -122.40379 - ...
## vertices: metric_graph_vertices
## VtEfirst: function ()
## Private:
## A: function (group = NULL, obs_to_vert = FALSE, drop_na = FALSE,
## add_vertices: function (PtE, tolerance = 1e-10, verbose)
## addinfo: FALSE
## clear_initial_info: function ()
## compute_degrees: function ()
## compute_laplacian_PtE: function (PtE, normalized = TRUE, verbose = verbose)
## compute_lengths: function (longlat, unit, crs, proj4string, which_longlat, vertex_unit,
## connected: TRUE
## coordinates_multiple_snaps: function (XY, tolerance, verbose = verbose, crs, proj4string,
## create_update_vertices: function ()
## crs: crs
## data: metric_graph_data, list
## edge_weights: tbl_df, tbl, data.frame
## find_edge_edge_points: function (tol, verbose, crs, proj4string, longlat, fact, which_longlat)
## find_mesh_bc: function ()
## get_edge_weights_internal: function (data.frame = FALSE)
## group_col: .group
## initial_edges_added: NULL
## initial_graph: metric_graph, R6
## kirchhoff_weights: NULL
## length_unit: km
## line_to_vertex: function (tolerance = 0, longlat = FALSE, fact, verbose, crs,
## longlat: TRUE
## merge_close_vertices: function (tolerance, fact)
## merge.all.deg2: function ()
## mesh_merge_deg2: function ()
## mesh_merge_outs: function ()
## move_V_first: function ()
## plot_2d: function (line_width = 0.1, marker_size = 1, vertex_color = "black",
## plot_3d: function (line_width = 1, marker_size = 1, vertex_color = "rgb(0,0,0)",
## proj4string: NULL
## prune_warning: FALSE
## pruned: FALSE
## PtE_to_mesh: function (PtE)
## ref_edges: 1 1552 2 5 3133 4 6 395 7 526 8 6466 9 10 11 12 38 13 14 ...
## remove_circles: function (threshold, verbose, longlat, unit, crs, proj4string,
## remove.first.deg2: function (res)
## set_first_weights: function (weights = rep(1, self$nE))
## set_petrov_matrices: function ()
## split_edge: function (Ei, t, tolerance = 0)
## temp_PtE: NULL
## tolerance: list
## transform: FALSE
## vertex_unit: degrees
## which_longlat: sf
## $ fem_mesh :List of 5
## ..$ c0:Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
## .. .. ..@ i : int [1:7169] 0 1 2 3 4 5 6 7 8 9 ...
## .. .. ..@ j : int [1:7169] 0 1 2 3 4 5 6 7 8 9 ...
## .. .. ..@ Dim : int [1:2] 7169 7169
## .. .. ..@ Dimnames:List of 2
## .. .. .. ..$ : NULL
## .. .. .. ..$ : NULL
## .. .. ..@ x : num [1:7169] 0.0797 0.1564 0.1117 0.0493 0.0614 ...
## .. .. ..@ factors : list()
## ..$ g1:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## .. .. ..@ i : int [1:24879] 0 1 985 0 1 1700 5207 5364 6858 7041 ...
## .. .. ..@ p : int [1:7170] 0 3 10 16 20 23 27 36 44 50 ...
## .. .. ..@ Dim : int [1:2] 7169 7169
## .. .. ..@ Dimnames:List of 2
## .. .. .. ..$ : NULL
## .. .. .. ..$ : NULL
## .. .. ..@ x : num [1:24879] 101.1 -55.1 -46.1 -55.1 216.4 ...
## .. .. ..@ factors : list()
## ..$ g2:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## .. .. ..@ i : int [1:56007] 0 1 985 986 1700 4722 5207 5364 6240 6858 ...
## .. .. ..@ p : int [1:7170] 0 11 22 34 43 52 61 76 89 99 ...
## .. .. ..@ Dim : int [1:2] 7169 7169
## .. .. ..@ Dimnames:List of 2
## .. .. .. ..$ : NULL
## .. .. .. ..$ : NULL
## .. .. ..@ x : num [1:56007] 172097 -145986 -232979 119006 15770 ...
## .. .. ..@ factors : list()
## ..$ g3:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## .. .. ..@ i : int [1:96755] 0 1 230 984 985 986 1700 1701 2280 2455 ...
## .. .. ..@ p : int [1:7170] 0 17 35 57 76 91 103 126 150 166 ...
## .. .. ..@ Dim : int [1:2] 7169 7169
## .. .. ..@ Dimnames:List of 2
## .. .. .. ..$ : NULL
## .. .. .. ..$ : NULL
## .. .. ..@ x : num [1:96755] 3.94e+08 -3.74e+08 -2.65e+07 -9.14e+07 -1.47e+09 ...
## .. .. ..@ factors : list()
## ..$ g4:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## .. .. ..@ i : int [1:147901] 0 1 230 983 984 985 986 989 1700 1701 ...
## .. .. ..@ p : int [1:7170] 0 28 52 89 122 151 171 209 249 274 ...
## .. .. ..@ Dim : int [1:2] 7169 7169
## .. .. ..@ Dimnames:List of 2
## .. .. .. ..$ : NULL
## .. .. .. ..$ : NULL
## .. .. ..@ x : num [1:147901] 1.41e+12 -1.06e+12 -2.25e+11 1.91e+11 -1.41e+12 ...
## .. .. ..@ factors : list()
## $ parameterization : chr "matern"
## $ n.spde : int 7169
## - attr(*, "class")= chr [1:3] "rspde_metric_graph" "inla_rspde" "inla.cgeneric"data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat,
repl = ".all",
loc_name = "loc")str(data_rspde_bru_stat)
## List of 4
## $ data :List of 8
## ..$ speed : num [1:14534] 0 12.9 24.1 32.2 0 ...
## ..$ SpeedLimit : num [1:14534] -0.101 -0.617 -0.617 -0.617 -0.927 ...
## ..$ .coord_x : num [1:14534] -122 -122 -122 -122 -122 ...
## ..$ .coord_y : num [1:14534] 37.8 37.8 37.8 37.8 37.8 ...
## ..$ .edge_number : num [1:14534] 1 4 6 6 9 14 14 14 18 20 ...
## ..$ .distance_on_edge: num [1:14534] 0.437 0.144 0.252 0.658 0.601 ...
## ..$ .group : chr [1:14534] "1" "1" "1" "1" ...
## ..$ loc : num [1:14534, 1:2] 1 4 6 6 9 14 14 14 18 20 ...
## ..- attr(*, "class")= chr [1:2] "metric_graph_data" "list"
## $ index:List of 3
## ..$ field : int [1:28676] 1 2 3 4 5 6 7 8 9 10 ...
## ..$ field.group: int [1:28676] 1 1 1 1 1 1 1 1 1 1 ...
## ..$ field.repl : int [1:28676] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "class")= chr [1:2] "inla_rspde_index" "list"
## ..- attr(*, "rspde.order")= num 0
## ..- attr(*, "integer_nu")= logi TRUE
## ..- attr(*, "n.mesh")= int 7169
## ..- attr(*, "name")= chr "field"
## ..- attr(*, "n.group")= int 1
## ..- attr(*, "n.repl")= int 4
## $ repl : chr [1:14534] "1" "1" "1" "1" ...
## $ basis:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## .. ..@ i : int [1:87204] 0 555 1905 1906 0 1050 1471 1472 1473 1905 ...
## .. ..@ p : int [1:86029] 0 4 11 11 11 11 12 13 15 15 ...
## .. ..@ Dim : int [1:2] 14534 86028
## .. ..@ Dimnames:List of 2
## .. .. ..$ : NULL
## .. .. ..$ : NULL
## .. ..@ x : num [1:87204] 0.563 0.0495 0.3595 0.7935 0.437 ...
## .. ..@ factors : list()cmp_stat = speed ~ -1 +
Intercept(1) +
SpeedLimit +
field(loc, model = rspde_model_stat,
replicate = data_rspde_bru_stat[["repl"]])
rspde_fit_stat <-
bru(cmp_stat,
data = data_rspde_bru_stat[["data"]],
family = "gaussian",
options = list(verbose = FALSE)
)
## Error in inla.inlaprogram.has.crashed() :
## The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
## If this does not help, please contact the developers at <help@r-inla.org>.
##
## *** inla.core.safe: inla.program has crashed: rerun to get better initial values. try=1/2
## Error in inla.inlaprogram.has.crashed() :
## The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
## If this does not help, please contact the developers at <help@r-inla.org>.
##
## *** inla.core.safe: inla.program has crashed: rerun to get better initial values. try=2/2
## Error in inla.inlaprogram.has.crashed() :
## The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
## If this does not help, please contact the developers at <help@r-inla.org>.str(rspde_fit_stat)
## Error in eval(expr, envir, enclos): object 'rspde_fit_stat' not foundstat.time.fin <- Sys.time()
print(stat.time.fin - stat.time.ini)## Time difference of 1.323035 minssummary(rspde_fit_stat)## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'rspde_fit_stat' not foundfit.rspde = rspde.result(rspde_fit_stat, "field", rspde_model_stat)## Error in eval(expr, envir, enclos): object 'rspde_fit_stat' not foundsummary(fit.rspde)## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'fit.rspde' not found1.2 Nonstationary model
- Observe that we are using the computed parameters from the stationary model as initial values for the nonstationary models.
nonstat.time.ini <- Sys.time()
################################################################################
############################# NON STATIONARY MODEL #############################
################################################################################
B.sigma = cbind(0, 1, 0, mesh$SpeedLimit, 0)
B.range = cbind(0, 0, 1, 0, mesh$SpeedLimit)
init.vec.theta = c(fit.rspde$summary.log.std.dev$mode,
fit.rspde$summary.log.range$mode,
rep(0, (ncol(B.sigma)-3)))
## Error in eval(expr, envir, enclos): object 'fit.rspde' not found
rspde_model_nonstat <- rspde.metric_graph(sf_graph,
start.theta = init.vec.theta,
theta.prior.mean = init.vec.theta,
B.sigma = B.sigma,
B.range = B.range,
parameterization = "matern",
nu.upper.bound = 1.5)
## Error in eval(expr, envir, enclos): object 'init.vec.theta' not foundstr(rspde_model_nonstat)
## Error in eval(expr, envir, enclos): object 'rspde_model_nonstat' not founddata_rspde_bru_nonstat <- graph_data_rspde(rspde_model_nonstat,
repl = ".all",
loc_name = "loc")
## Error in eval(expr, envir, enclos): object 'rspde_model_nonstat' not foundstr(data_rspde_bru_nonstat)
## Error in eval(expr, envir, enclos): object 'data_rspde_bru_nonstat' not foundcmp_nonstat = speed ~ -1 +
Intercept(1) +
SpeedLimit +
field(loc, model = rspde_model_nonstat,
replicate = data_rspde_bru_nonstat[["repl"]])
rspde_fit_nonstat <-
bru(cmp_nonstat,
data = data_rspde_bru_nonstat[["data"]],
family = "gaussian",
options = list(verbose = FALSE)
)
## Error in eval(expr, envir, enclos): object 'data_rspde_bru_nonstat' not foundstr(rspde_fit_nonstat)
## Error in eval(expr, envir, enclos): object 'rspde_fit_nonstat' not foundnonstat.time.fin <- Sys.time()
print(nonstat.time.fin - nonstat.time.ini)## Time difference of 0.02544212 secssummary(rspde_fit_nonstat)## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'rspde_fit_nonstat' not foundsummary(rspde.result(rspde_fit_nonstat, "field", rspde_model_nonstat))## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'rspde_model_nonstat' not found1.3 Crossvalidation 1
#load(here("Models_output/distmatrixfixed.RData"))
points = data %>%
as.data.frame() %>%
st_as_sf(coords = c(".coord_x", ".coord_y"), crs = 4326) %>%
mutate(., index = 1:nrow(.)) %>%
st_drop_geometry() %>%
dplyr:::select(speed, .group, index) %>%
mutate(.group = as.numeric(.group)) %>%
group_by(.group) %>%
mutate(indexingroup = seq_len(n())) %>%
ungroup()
distance = seq(from = 0, to = 400, by = 20)/1000
#save(GROUPS, file = here("Models_output/GROUPS_for_window_case.RData"))load(here("Models_output/GROUPS_for_window_case.RData"))
mse.stat <- mse.nonstat <- ls.stat <- ls.nonstat <- rep(0,length(distance))
# cross-validation for-loop
for (j in 1:length(distance)) {
print(j)
# cross-validation of the stationary model
cv.stat <- inla.group.cv(rspde_fit_stat, groups = GROUPS[[j]])
# cross-validation of the nonstationary model
cv.nonstat <- inla.group.cv(rspde_fit_nonstat, groups = GROUPS[[j]])
# obtain MSE and LS
mse.stat[j] <- mean((cv.stat$mean - points$speed)^2)
mse.nonstat[j] <- mean((cv.nonstat$mean - points$speed)^2)
ls.stat[j] <- mean(log(cv.stat$cv))
ls.nonstat[j] <- mean(log(cv.nonstat$cv))
}
## [1] 1
## Error in eval(expr, envir, enclos): object 'rspde_fit_stat' not found
## plot results
par(mfrow = c(2,2), family = "Palatino")
# Plot MSE
plot(distance, mse.stat, main = "MSE", ylim = c(min(mse.nonstat, mse.stat), max(mse.nonstat, mse.stat)),
type = "l", ylab = "MSE", xlab = "distance in m", col = "black")
lines(distance, mse.nonstat, col = "blue")
legend("bottomright", legend = c("Stationary", "Non-stationary"), col = c("black", "blue"), lty = 1)
# Plot log-score
plot(distance, -ls.stat, main = "log-score", ylim = c(min(-ls.nonstat, -ls.stat), max(-ls.nonstat, -ls.stat)),
type = "l", ylab = "log-score", xlab = "distance in m", col = "black")
lines(distance, -ls.nonstat, col = "blue")
legend("bottomright", legend = c("Stationary", "Non-stationary"), col = c("black", "blue"), lty = 1)
save.image(here(paste0("Models_output/", rmarkdown::metadata$title, ".RData")))